Linear Regression: Pollution Example

Before beginning the analysis, it is helpful to start by visualizing the data.

Id maxO3 T9 T12 T15 Ne9 Ne12 Ne15 Vx9 Vx12 Vx15 maxO3v wind rain
20010601 87 15.6 18.5 18.4 4 4 8 0.69 -1.71 -0.69 84 Nord Sec
20010602 82 17.0 18.4 17.7 5 5 7 -4.33 -4.00 -3.00 87 Nord Sec
20010603 92 15.3 17.6 19.5 2 5 4 2.95 1.88 0.52 82 Est Sec
20010604 114 16.2 19.7 22.5 1 1 0 0.98 0.35 -0.17 92 Nord Sec
20010605 94 17.4 20.5 20.4 8 8 7 -0.50 -2.95 -4.33 114 Ouest Sec
20010606 80 17.7 19.8 18.3 6 6 7 -5.64 -5.00 -6.00 94 Ouest Pluie

Correlation

## corrplot 0.95 loaded

library("scatterplot3d")
scatterplot3d(ozone[,"T12"],ozone[,"Vx12"],ozone[,"maxO3"],type="h", pch=16,box=FALSE,xlab="Température à 12h",ylab="Vent",zlab="Ozone")

## When the points are high, this means that the rate of ozone is higher in the air. We can notice that the ozone increases as the temperature is increases also.

Data

##         Id maxO3   T9  T12  T15 Ne9 Ne12 Ne15     Vx9    Vx12    Vx15 maxO3v
## 1 20010601    87 15.6 18.5 18.4   4    4    8  0.6946 -1.7101 -0.6946     84
## 2 20010602    82 17.0 18.4 17.7   5    5    7 -4.3301 -4.0000 -3.0000     87
## 3 20010603    92 15.3 17.6 19.5   2    5    4  2.9544  1.8794  0.5209     82
## 4 20010604   114 16.2 19.7 22.5   1    1    0  0.9848  0.3473 -0.1736     92
## 5 20010605    94 17.4 20.5 20.4   8    8    7 -0.5000 -2.9544 -4.3301    114
## 6 20010606    80 17.7 19.8 18.3   6    6    7 -5.6382 -5.0000 -6.0000     94
##    wind  rain
## 1  Nord   Sec
## 2  Nord   Sec
## 3   Est   Sec
## 4  Nord   Sec
## 5 Ouest   Sec
## 6 Ouest Pluie

PCA

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
##   Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Temperature and wind components contribute in 56.45% in the first dimension, while the cloud cover contribute at 17.16% in the second dimension.

When the wind blows in the west direction, there is a cloud cover.

Example of pollution

We consider a linear model between ozone, wind, and temperature.

modele1 <- lm(maxO3~Vx12+T12,data=ozone)
summary(modele1)
## 
## Call:
## lm(formula = maxO3 ~ Vx12 + T12, data = ozone)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -33.08 -11.99  -1.20  10.67  45.67 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -14.4242     9.3943  -1.535  0.12758    
## Vx12          2.0742     0.5987   3.465  0.00076 ***
## T12           5.0202     0.4140  12.125  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.75 on 109 degrees of freedom
## Multiple R-squared:  0.6533, Adjusted R-squared:  0.6469 
## F-statistic: 102.7 on 2 and 109 DF,  p-value: < 2.2e-16
anova(modele1)
## Analysis of Variance Table
## 
## Response: maxO3
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Vx12        1  16367   16367  58.338 9.074e-12 ***
## T12         1  41244   41244 147.010 < 2.2e-16 ***
## Residuals 109  30580     281                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Only the wind component and temperature at 12PM are significant. The intercept is not significant.The model only explain about 64.7% of the variance which is not good enough.

Model without intercept

## 
## Call:
## lm(formula = maxO3 ~ -1 + Vx12 + T12, data = ozone)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -32.135 -12.914  -1.394  10.020  48.525 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## Vx12   2.4412     0.5523    4.42 2.32e-05 ***
## T12    4.3967     0.0811   54.22  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.85 on 110 degrees of freedom
## Multiple R-squared:  0.9688, Adjusted R-squared:  0.9682 
## F-statistic:  1708 on 2 and 110 DF,  p-value: < 2.2e-16

When we remove the intercept which was previously not significant, the new model explain about 96% of the variance which is good.

Without wind

## 
## Call:
## lm(formula = maxO3 ~ Ne12, data = ozone)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.020 -14.487  -5.115  12.571  66.470 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 130.0201     4.9807  26.105  < 2e-16 ***
## Ne12         -7.9150     0.9042  -8.753 2.77e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.74 on 110 degrees of freedom
## Multiple R-squared:  0.4106, Adjusted R-squared:  0.4052 
## F-statistic: 76.62 on 1 and 110 DF,  p-value: 2.769e-14

The adjusted R-squared explains only 40% of the variance which is very low.

## 
## Call:
## lm(formula = maxO3 ~ Vx9 + Vx12 + Vx15 + Ne9 + Ne12 + Ne15 + 
##     T9 + T12 + T15 + maxO3v, data = ozone)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -53.566  -8.727  -0.403   7.599  39.458 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.24442   13.47190   0.909   0.3656    
## Vx9          0.94791    0.91228   1.039   0.3013    
## Vx12         0.03120    1.05523   0.030   0.9765    
## Vx15         0.41859    0.91568   0.457   0.6486    
## Ne9         -2.18909    0.93824  -2.333   0.0216 *  
## Ne12        -0.42102    1.36766  -0.308   0.7588    
## Ne15         0.18373    1.00279   0.183   0.8550    
## T9          -0.01901    1.12515  -0.017   0.9866    
## T12          2.22115    1.43294   1.550   0.1243    
## T15          0.55853    1.14464   0.488   0.6266    
## maxO3v       0.35198    0.06289   5.597 1.88e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.36 on 101 degrees of freedom
## Multiple R-squared:  0.7638, Adjusted R-squared:  0.7405 
## F-statistic: 32.67 on 10 and 101 DF,  p-value: < 2.2e-16

Only the cloud cover(Ne9) at 9AM and previous max ozone(max03v) are significant and the model explains 74% of the variance. Even though this is a good adjusted squared, the model is very complex as it has a lot of variables, most of which are not significant.

The most important variable is the temperature at 12PM(T12) and the previous max ozone(maxo3v). And the cloud cover at 9AM(Ne9) in the third position.

R2 Choice

Keep Vx9, Ne9, T12, maxO3v

plot(choix,scale="adjr2")

By BIC

Keep Intercept Vx9, Ne9, T12, maxO3v

plot(choix,scale="bic")

If p>>n?

Let us consider the data on the contents of fat, sugar, flour, and water in biscuits (Osbone et al., 1984) for 72 biscuits, available in the “ppls” package. For these 72 biscuits, the compositions in fat, sugar, flour, and water are measured using a classical approach, while the spectrum is observed over all wavelengths between 1100 and 2498 nanometers, regularly spaced every 2 nanometers. We therefore have 700 observed values, or potentially explanatory variables, for each biscuit dough sample.

Typically, this study is conducted in a very high-dimensional context with \(p >> n\). On these uncooked biscuits, absorbance is measured for wavelengths between 1100 and 2498 nanometers, regularly spaced every 2 nanometers. There are thus 700 potentially explanatory variables.

Let us consider the percentage of sugars. We have \(p = 700\) variables for \(n = 40\) individuals. The classical least squares estimator is not defined.

The classical least squares estimator is not defined.

Lectrure des données

library(ppls)
data(cookie)
# extraire le taux de sucrose et les spectres
cook = data.frame(cookie$constituents$sucrose,cookie$NIR)
names(cook)= c("sucrose",paste("X",1:700,sep=""))

Data

# Response
hist(cook[,"sucrose"])

## The sucrose percentage is between 8% and 22%. The distribution is most concentrated between 10% and 20%. This shows moderate variability in sugar content across the biscuits samples

# spectres by sucrose intensity
coul=rainbow(20)[as.integer(as.factor
    (as.integer(cook[,1]-10)))]
ts.plot(t(cook[,-1]),col=coul) 

## The spectral plot shows us many colors mean that we have to have to use multivariate method instead of univariate method.

MCO estimation not valide

cook.app=cook[1:40,]
summary(lm(sucrose ~ ., data=cook.app))
## 
## Call:
## lm(formula = sucrose ~ ., data = cook.app)
## 
## Residuals:
## ALL 40 residuals are 0: no residual degrees of freedom!
## 
## Coefficients: (661 not defined because of singularities)
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)     12.17        NaN     NaN      NaN
## X1           13724.43        NaN     NaN      NaN
## X2          -31660.99        NaN     NaN      NaN
## X3           25153.86        NaN     NaN      NaN
## X4          -13071.94        NaN     NaN      NaN
## X5             244.71        NaN     NaN      NaN
## X6           24773.95        NaN     NaN      NaN
## X7          -26093.74        NaN     NaN      NaN
## X8           12801.63        NaN     NaN      NaN
## X9           -4709.51        NaN     NaN      NaN
## X10          -4891.93        NaN     NaN      NaN
## X11          31636.00        NaN     NaN      NaN
## X12         -21089.99        NaN     NaN      NaN
## X13          -1917.31        NaN     NaN      NaN
## X14         -14403.37        NaN     NaN      NaN
## X15          15061.66        NaN     NaN      NaN
## X16          22386.40        NaN     NaN      NaN
## X17         -54867.14        NaN     NaN      NaN
## X18            287.33        NaN     NaN      NaN
## X19          -8562.09        NaN     NaN      NaN
## X20          34838.19        NaN     NaN      NaN
## X21          17488.40        NaN     NaN      NaN
## X22         -23338.41        NaN     NaN      NaN
## X23          20792.72        NaN     NaN      NaN
## X24         -32664.20        NaN     NaN      NaN
## X25          35744.76        NaN     NaN      NaN
## X26          12087.43        NaN     NaN      NaN
## X27         -13153.13        NaN     NaN      NaN
## X28         -28838.43        NaN     NaN      NaN
## X29          22346.68        NaN     NaN      NaN
## X30         -17483.16        NaN     NaN      NaN
## X31          28184.57        NaN     NaN      NaN
## X32         -25085.04        NaN     NaN      NaN
## X33          -6570.85        NaN     NaN      NaN
## X34          -8131.31        NaN     NaN      NaN
## X35          14008.36        NaN     NaN      NaN
## X36          -6829.77        NaN     NaN      NaN
## X37          10198.76        NaN     NaN      NaN
## X38         -13024.01        NaN     NaN      NaN
## X39          14510.81        NaN     NaN      NaN
## X40                NA         NA      NA       NA
## X41                NA         NA      NA       NA
## X42                NA         NA      NA       NA
## X43                NA         NA      NA       NA
## X44                NA         NA      NA       NA
## X45                NA         NA      NA       NA
## X46                NA         NA      NA       NA
## X47                NA         NA      NA       NA
## X48                NA         NA      NA       NA
## X49                NA         NA      NA       NA
## X50                NA         NA      NA       NA
## X51                NA         NA      NA       NA
## X52                NA         NA      NA       NA
## X53                NA         NA      NA       NA
## X54                NA         NA      NA       NA
## X55                NA         NA      NA       NA
## X56                NA         NA      NA       NA
## X57                NA         NA      NA       NA
## X58                NA         NA      NA       NA
## X59                NA         NA      NA       NA
## X60                NA         NA      NA       NA
## X61                NA         NA      NA       NA
## X62                NA         NA      NA       NA
## X63                NA         NA      NA       NA
## X64                NA         NA      NA       NA
## X65                NA         NA      NA       NA
## X66                NA         NA      NA       NA
## X67                NA         NA      NA       NA
## X68                NA         NA      NA       NA
## X69                NA         NA      NA       NA
## X70                NA         NA      NA       NA
## X71                NA         NA      NA       NA
## X72                NA         NA      NA       NA
## X73                NA         NA      NA       NA
## X74                NA         NA      NA       NA
## X75                NA         NA      NA       NA
## X76                NA         NA      NA       NA
## X77                NA         NA      NA       NA
## X78                NA         NA      NA       NA
## X79                NA         NA      NA       NA
## X80                NA         NA      NA       NA
## X81                NA         NA      NA       NA
## X82                NA         NA      NA       NA
## X83                NA         NA      NA       NA
## X84                NA         NA      NA       NA
## X85                NA         NA      NA       NA
## X86                NA         NA      NA       NA
## X87                NA         NA      NA       NA
## X88                NA         NA      NA       NA
## X89                NA         NA      NA       NA
## X90                NA         NA      NA       NA
## X91                NA         NA      NA       NA
## X92                NA         NA      NA       NA
## X93                NA         NA      NA       NA
## X94                NA         NA      NA       NA
## X95                NA         NA      NA       NA
## X96                NA         NA      NA       NA
## X97                NA         NA      NA       NA
## X98                NA         NA      NA       NA
## X99                NA         NA      NA       NA
## X100               NA         NA      NA       NA
## X101               NA         NA      NA       NA
## X102               NA         NA      NA       NA
## X103               NA         NA      NA       NA
## X104               NA         NA      NA       NA
## X105               NA         NA      NA       NA
## X106               NA         NA      NA       NA
## X107               NA         NA      NA       NA
## X108               NA         NA      NA       NA
## X109               NA         NA      NA       NA
## X110               NA         NA      NA       NA
## X111               NA         NA      NA       NA
## X112               NA         NA      NA       NA
## X113               NA         NA      NA       NA
## X114               NA         NA      NA       NA
## X115               NA         NA      NA       NA
## X116               NA         NA      NA       NA
## X117               NA         NA      NA       NA
## X118               NA         NA      NA       NA
## X119               NA         NA      NA       NA
## X120               NA         NA      NA       NA
## X121               NA         NA      NA       NA
## X122               NA         NA      NA       NA
## X123               NA         NA      NA       NA
## X124               NA         NA      NA       NA
## X125               NA         NA      NA       NA
## X126               NA         NA      NA       NA
## X127               NA         NA      NA       NA
## X128               NA         NA      NA       NA
## X129               NA         NA      NA       NA
## X130               NA         NA      NA       NA
## X131               NA         NA      NA       NA
## X132               NA         NA      NA       NA
## X133               NA         NA      NA       NA
## X134               NA         NA      NA       NA
## X135               NA         NA      NA       NA
## X136               NA         NA      NA       NA
## X137               NA         NA      NA       NA
## X138               NA         NA      NA       NA
## X139               NA         NA      NA       NA
## X140               NA         NA      NA       NA
## X141               NA         NA      NA       NA
## X142               NA         NA      NA       NA
## X143               NA         NA      NA       NA
## X144               NA         NA      NA       NA
## X145               NA         NA      NA       NA
## X146               NA         NA      NA       NA
## X147               NA         NA      NA       NA
## X148               NA         NA      NA       NA
## X149               NA         NA      NA       NA
## X150               NA         NA      NA       NA
## X151               NA         NA      NA       NA
## X152               NA         NA      NA       NA
## X153               NA         NA      NA       NA
## X154               NA         NA      NA       NA
## X155               NA         NA      NA       NA
## X156               NA         NA      NA       NA
## X157               NA         NA      NA       NA
## X158               NA         NA      NA       NA
## X159               NA         NA      NA       NA
## X160               NA         NA      NA       NA
## X161               NA         NA      NA       NA
## X162               NA         NA      NA       NA
## X163               NA         NA      NA       NA
## X164               NA         NA      NA       NA
## X165               NA         NA      NA       NA
## X166               NA         NA      NA       NA
## X167               NA         NA      NA       NA
## X168               NA         NA      NA       NA
## X169               NA         NA      NA       NA
## X170               NA         NA      NA       NA
## X171               NA         NA      NA       NA
## X172               NA         NA      NA       NA
## X173               NA         NA      NA       NA
## X174               NA         NA      NA       NA
## X175               NA         NA      NA       NA
## X176               NA         NA      NA       NA
## X177               NA         NA      NA       NA
## X178               NA         NA      NA       NA
## X179               NA         NA      NA       NA
## X180               NA         NA      NA       NA
## X181               NA         NA      NA       NA
## X182               NA         NA      NA       NA
## X183               NA         NA      NA       NA
## X184               NA         NA      NA       NA
## X185               NA         NA      NA       NA
## X186               NA         NA      NA       NA
## X187               NA         NA      NA       NA
## X188               NA         NA      NA       NA
## X189               NA         NA      NA       NA
## X190               NA         NA      NA       NA
## X191               NA         NA      NA       NA
## X192               NA         NA      NA       NA
## X193               NA         NA      NA       NA
## X194               NA         NA      NA       NA
## X195               NA         NA      NA       NA
## X196               NA         NA      NA       NA
## X197               NA         NA      NA       NA
## X198               NA         NA      NA       NA
## X199               NA         NA      NA       NA
## X200               NA         NA      NA       NA
## X201               NA         NA      NA       NA
## X202               NA         NA      NA       NA
## X203               NA         NA      NA       NA
## X204               NA         NA      NA       NA
## X205               NA         NA      NA       NA
## X206               NA         NA      NA       NA
## X207               NA         NA      NA       NA
## X208               NA         NA      NA       NA
## X209               NA         NA      NA       NA
## X210               NA         NA      NA       NA
## X211               NA         NA      NA       NA
## X212               NA         NA      NA       NA
## X213               NA         NA      NA       NA
## X214               NA         NA      NA       NA
## X215               NA         NA      NA       NA
## X216               NA         NA      NA       NA
## X217               NA         NA      NA       NA
## X218               NA         NA      NA       NA
## X219               NA         NA      NA       NA
## X220               NA         NA      NA       NA
## X221               NA         NA      NA       NA
## X222               NA         NA      NA       NA
## X223               NA         NA      NA       NA
## X224               NA         NA      NA       NA
## X225               NA         NA      NA       NA
## X226               NA         NA      NA       NA
## X227               NA         NA      NA       NA
## X228               NA         NA      NA       NA
## X229               NA         NA      NA       NA
## X230               NA         NA      NA       NA
## X231               NA         NA      NA       NA
## X232               NA         NA      NA       NA
## X233               NA         NA      NA       NA
## X234               NA         NA      NA       NA
## X235               NA         NA      NA       NA
## X236               NA         NA      NA       NA
## X237               NA         NA      NA       NA
## X238               NA         NA      NA       NA
## X239               NA         NA      NA       NA
## X240               NA         NA      NA       NA
## X241               NA         NA      NA       NA
## X242               NA         NA      NA       NA
## X243               NA         NA      NA       NA
## X244               NA         NA      NA       NA
## X245               NA         NA      NA       NA
## X246               NA         NA      NA       NA
## X247               NA         NA      NA       NA
## X248               NA         NA      NA       NA
## X249               NA         NA      NA       NA
## X250               NA         NA      NA       NA
## X251               NA         NA      NA       NA
## X252               NA         NA      NA       NA
## X253               NA         NA      NA       NA
## X254               NA         NA      NA       NA
## X255               NA         NA      NA       NA
## X256               NA         NA      NA       NA
## X257               NA         NA      NA       NA
## X258               NA         NA      NA       NA
## X259               NA         NA      NA       NA
## X260               NA         NA      NA       NA
## X261               NA         NA      NA       NA
## X262               NA         NA      NA       NA
## X263               NA         NA      NA       NA
## X264               NA         NA      NA       NA
## X265               NA         NA      NA       NA
## X266               NA         NA      NA       NA
## X267               NA         NA      NA       NA
## X268               NA         NA      NA       NA
## X269               NA         NA      NA       NA
## X270               NA         NA      NA       NA
## X271               NA         NA      NA       NA
## X272               NA         NA      NA       NA
## X273               NA         NA      NA       NA
## X274               NA         NA      NA       NA
## X275               NA         NA      NA       NA
## X276               NA         NA      NA       NA
## X277               NA         NA      NA       NA
## X278               NA         NA      NA       NA
## X279               NA         NA      NA       NA
## X280               NA         NA      NA       NA
## X281               NA         NA      NA       NA
## X282               NA         NA      NA       NA
## X283               NA         NA      NA       NA
## X284               NA         NA      NA       NA
## X285               NA         NA      NA       NA
## X286               NA         NA      NA       NA
## X287               NA         NA      NA       NA
## X288               NA         NA      NA       NA
## X289               NA         NA      NA       NA
## X290               NA         NA      NA       NA
## X291               NA         NA      NA       NA
## X292               NA         NA      NA       NA
## X293               NA         NA      NA       NA
## X294               NA         NA      NA       NA
## X295               NA         NA      NA       NA
## X296               NA         NA      NA       NA
## X297               NA         NA      NA       NA
## X298               NA         NA      NA       NA
## X299               NA         NA      NA       NA
## X300               NA         NA      NA       NA
## X301               NA         NA      NA       NA
## X302               NA         NA      NA       NA
## X303               NA         NA      NA       NA
## X304               NA         NA      NA       NA
## X305               NA         NA      NA       NA
## X306               NA         NA      NA       NA
## X307               NA         NA      NA       NA
## X308               NA         NA      NA       NA
## X309               NA         NA      NA       NA
## X310               NA         NA      NA       NA
## X311               NA         NA      NA       NA
## X312               NA         NA      NA       NA
## X313               NA         NA      NA       NA
## X314               NA         NA      NA       NA
## X315               NA         NA      NA       NA
## X316               NA         NA      NA       NA
## X317               NA         NA      NA       NA
## X318               NA         NA      NA       NA
## X319               NA         NA      NA       NA
## X320               NA         NA      NA       NA
## X321               NA         NA      NA       NA
## X322               NA         NA      NA       NA
## X323               NA         NA      NA       NA
## X324               NA         NA      NA       NA
## X325               NA         NA      NA       NA
## X326               NA         NA      NA       NA
## X327               NA         NA      NA       NA
## X328               NA         NA      NA       NA
## X329               NA         NA      NA       NA
## X330               NA         NA      NA       NA
## X331               NA         NA      NA       NA
## X332               NA         NA      NA       NA
## X333               NA         NA      NA       NA
## X334               NA         NA      NA       NA
## X335               NA         NA      NA       NA
## X336               NA         NA      NA       NA
## X337               NA         NA      NA       NA
## X338               NA         NA      NA       NA
## X339               NA         NA      NA       NA
## X340               NA         NA      NA       NA
## X341               NA         NA      NA       NA
## X342               NA         NA      NA       NA
## X343               NA         NA      NA       NA
## X344               NA         NA      NA       NA
## X345               NA         NA      NA       NA
## X346               NA         NA      NA       NA
## X347               NA         NA      NA       NA
## X348               NA         NA      NA       NA
## X349               NA         NA      NA       NA
## X350               NA         NA      NA       NA
## X351               NA         NA      NA       NA
## X352               NA         NA      NA       NA
## X353               NA         NA      NA       NA
## X354               NA         NA      NA       NA
## X355               NA         NA      NA       NA
## X356               NA         NA      NA       NA
## X357               NA         NA      NA       NA
## X358               NA         NA      NA       NA
## X359               NA         NA      NA       NA
## X360               NA         NA      NA       NA
## X361               NA         NA      NA       NA
## X362               NA         NA      NA       NA
## X363               NA         NA      NA       NA
## X364               NA         NA      NA       NA
## X365               NA         NA      NA       NA
## X366               NA         NA      NA       NA
## X367               NA         NA      NA       NA
## X368               NA         NA      NA       NA
## X369               NA         NA      NA       NA
## X370               NA         NA      NA       NA
## X371               NA         NA      NA       NA
## X372               NA         NA      NA       NA
## X373               NA         NA      NA       NA
## X374               NA         NA      NA       NA
## X375               NA         NA      NA       NA
## X376               NA         NA      NA       NA
## X377               NA         NA      NA       NA
## X378               NA         NA      NA       NA
## X379               NA         NA      NA       NA
## X380               NA         NA      NA       NA
## X381               NA         NA      NA       NA
## X382               NA         NA      NA       NA
## X383               NA         NA      NA       NA
## X384               NA         NA      NA       NA
## X385               NA         NA      NA       NA
## X386               NA         NA      NA       NA
## X387               NA         NA      NA       NA
## X388               NA         NA      NA       NA
## X389               NA         NA      NA       NA
## X390               NA         NA      NA       NA
## X391               NA         NA      NA       NA
## X392               NA         NA      NA       NA
## X393               NA         NA      NA       NA
## X394               NA         NA      NA       NA
## X395               NA         NA      NA       NA
## X396               NA         NA      NA       NA
## X397               NA         NA      NA       NA
## X398               NA         NA      NA       NA
## X399               NA         NA      NA       NA
## X400               NA         NA      NA       NA
## X401               NA         NA      NA       NA
## X402               NA         NA      NA       NA
## X403               NA         NA      NA       NA
## X404               NA         NA      NA       NA
## X405               NA         NA      NA       NA
## X406               NA         NA      NA       NA
## X407               NA         NA      NA       NA
## X408               NA         NA      NA       NA
## X409               NA         NA      NA       NA
## X410               NA         NA      NA       NA
## X411               NA         NA      NA       NA
## X412               NA         NA      NA       NA
## X413               NA         NA      NA       NA
## X414               NA         NA      NA       NA
## X415               NA         NA      NA       NA
## X416               NA         NA      NA       NA
## X417               NA         NA      NA       NA
## X418               NA         NA      NA       NA
## X419               NA         NA      NA       NA
## X420               NA         NA      NA       NA
## X421               NA         NA      NA       NA
## X422               NA         NA      NA       NA
## X423               NA         NA      NA       NA
## X424               NA         NA      NA       NA
## X425               NA         NA      NA       NA
## X426               NA         NA      NA       NA
## X427               NA         NA      NA       NA
## X428               NA         NA      NA       NA
## X429               NA         NA      NA       NA
## X430               NA         NA      NA       NA
## X431               NA         NA      NA       NA
## X432               NA         NA      NA       NA
## X433               NA         NA      NA       NA
## X434               NA         NA      NA       NA
## X435               NA         NA      NA       NA
## X436               NA         NA      NA       NA
## X437               NA         NA      NA       NA
## X438               NA         NA      NA       NA
## X439               NA         NA      NA       NA
## X440               NA         NA      NA       NA
## X441               NA         NA      NA       NA
## X442               NA         NA      NA       NA
## X443               NA         NA      NA       NA
## X444               NA         NA      NA       NA
## X445               NA         NA      NA       NA
## X446               NA         NA      NA       NA
## X447               NA         NA      NA       NA
## X448               NA         NA      NA       NA
## X449               NA         NA      NA       NA
## X450               NA         NA      NA       NA
## X451               NA         NA      NA       NA
## X452               NA         NA      NA       NA
## X453               NA         NA      NA       NA
## X454               NA         NA      NA       NA
## X455               NA         NA      NA       NA
## X456               NA         NA      NA       NA
## X457               NA         NA      NA       NA
## X458               NA         NA      NA       NA
## X459               NA         NA      NA       NA
## X460               NA         NA      NA       NA
## X461               NA         NA      NA       NA
## X462               NA         NA      NA       NA
## X463               NA         NA      NA       NA
## X464               NA         NA      NA       NA
## X465               NA         NA      NA       NA
## X466               NA         NA      NA       NA
## X467               NA         NA      NA       NA
## X468               NA         NA      NA       NA
## X469               NA         NA      NA       NA
## X470               NA         NA      NA       NA
## X471               NA         NA      NA       NA
## X472               NA         NA      NA       NA
## X473               NA         NA      NA       NA
## X474               NA         NA      NA       NA
## X475               NA         NA      NA       NA
## X476               NA         NA      NA       NA
## X477               NA         NA      NA       NA
## X478               NA         NA      NA       NA
## X479               NA         NA      NA       NA
## X480               NA         NA      NA       NA
## X481               NA         NA      NA       NA
## X482               NA         NA      NA       NA
## X483               NA         NA      NA       NA
## X484               NA         NA      NA       NA
## X485               NA         NA      NA       NA
## X486               NA         NA      NA       NA
## X487               NA         NA      NA       NA
## X488               NA         NA      NA       NA
## X489               NA         NA      NA       NA
## X490               NA         NA      NA       NA
## X491               NA         NA      NA       NA
## X492               NA         NA      NA       NA
## X493               NA         NA      NA       NA
## X494               NA         NA      NA       NA
## X495               NA         NA      NA       NA
## X496               NA         NA      NA       NA
## X497               NA         NA      NA       NA
## X498               NA         NA      NA       NA
## X499               NA         NA      NA       NA
## X500               NA         NA      NA       NA
## X501               NA         NA      NA       NA
## X502               NA         NA      NA       NA
## X503               NA         NA      NA       NA
## X504               NA         NA      NA       NA
## X505               NA         NA      NA       NA
## X506               NA         NA      NA       NA
## X507               NA         NA      NA       NA
## X508               NA         NA      NA       NA
## X509               NA         NA      NA       NA
## X510               NA         NA      NA       NA
## X511               NA         NA      NA       NA
## X512               NA         NA      NA       NA
## X513               NA         NA      NA       NA
## X514               NA         NA      NA       NA
## X515               NA         NA      NA       NA
## X516               NA         NA      NA       NA
## X517               NA         NA      NA       NA
## X518               NA         NA      NA       NA
## X519               NA         NA      NA       NA
## X520               NA         NA      NA       NA
## X521               NA         NA      NA       NA
## X522               NA         NA      NA       NA
## X523               NA         NA      NA       NA
## X524               NA         NA      NA       NA
## X525               NA         NA      NA       NA
## X526               NA         NA      NA       NA
## X527               NA         NA      NA       NA
## X528               NA         NA      NA       NA
## X529               NA         NA      NA       NA
## X530               NA         NA      NA       NA
## X531               NA         NA      NA       NA
## X532               NA         NA      NA       NA
## X533               NA         NA      NA       NA
## X534               NA         NA      NA       NA
## X535               NA         NA      NA       NA
## X536               NA         NA      NA       NA
## X537               NA         NA      NA       NA
## X538               NA         NA      NA       NA
## X539               NA         NA      NA       NA
## X540               NA         NA      NA       NA
## X541               NA         NA      NA       NA
## X542               NA         NA      NA       NA
## X543               NA         NA      NA       NA
## X544               NA         NA      NA       NA
## X545               NA         NA      NA       NA
## X546               NA         NA      NA       NA
## X547               NA         NA      NA       NA
## X548               NA         NA      NA       NA
## X549               NA         NA      NA       NA
## X550               NA         NA      NA       NA
## X551               NA         NA      NA       NA
## X552               NA         NA      NA       NA
## X553               NA         NA      NA       NA
## X554               NA         NA      NA       NA
## X555               NA         NA      NA       NA
## X556               NA         NA      NA       NA
## X557               NA         NA      NA       NA
## X558               NA         NA      NA       NA
## X559               NA         NA      NA       NA
## X560               NA         NA      NA       NA
## X561               NA         NA      NA       NA
## X562               NA         NA      NA       NA
## X563               NA         NA      NA       NA
## X564               NA         NA      NA       NA
## X565               NA         NA      NA       NA
## X566               NA         NA      NA       NA
## X567               NA         NA      NA       NA
## X568               NA         NA      NA       NA
## X569               NA         NA      NA       NA
## X570               NA         NA      NA       NA
## X571               NA         NA      NA       NA
## X572               NA         NA      NA       NA
## X573               NA         NA      NA       NA
## X574               NA         NA      NA       NA
## X575               NA         NA      NA       NA
## X576               NA         NA      NA       NA
## X577               NA         NA      NA       NA
## X578               NA         NA      NA       NA
## X579               NA         NA      NA       NA
## X580               NA         NA      NA       NA
## X581               NA         NA      NA       NA
## X582               NA         NA      NA       NA
## X583               NA         NA      NA       NA
## X584               NA         NA      NA       NA
## X585               NA         NA      NA       NA
## X586               NA         NA      NA       NA
## X587               NA         NA      NA       NA
## X588               NA         NA      NA       NA
## X589               NA         NA      NA       NA
## X590               NA         NA      NA       NA
## X591               NA         NA      NA       NA
## X592               NA         NA      NA       NA
## X593               NA         NA      NA       NA
## X594               NA         NA      NA       NA
## X595               NA         NA      NA       NA
## X596               NA         NA      NA       NA
## X597               NA         NA      NA       NA
## X598               NA         NA      NA       NA
## X599               NA         NA      NA       NA
## X600               NA         NA      NA       NA
## X601               NA         NA      NA       NA
## X602               NA         NA      NA       NA
## X603               NA         NA      NA       NA
## X604               NA         NA      NA       NA
## X605               NA         NA      NA       NA
## X606               NA         NA      NA       NA
## X607               NA         NA      NA       NA
## X608               NA         NA      NA       NA
## X609               NA         NA      NA       NA
## X610               NA         NA      NA       NA
## X611               NA         NA      NA       NA
## X612               NA         NA      NA       NA
## X613               NA         NA      NA       NA
## X614               NA         NA      NA       NA
## X615               NA         NA      NA       NA
## X616               NA         NA      NA       NA
## X617               NA         NA      NA       NA
## X618               NA         NA      NA       NA
## X619               NA         NA      NA       NA
## X620               NA         NA      NA       NA
## X621               NA         NA      NA       NA
## X622               NA         NA      NA       NA
## X623               NA         NA      NA       NA
## X624               NA         NA      NA       NA
## X625               NA         NA      NA       NA
## X626               NA         NA      NA       NA
## X627               NA         NA      NA       NA
## X628               NA         NA      NA       NA
## X629               NA         NA      NA       NA
## X630               NA         NA      NA       NA
## X631               NA         NA      NA       NA
## X632               NA         NA      NA       NA
## X633               NA         NA      NA       NA
## X634               NA         NA      NA       NA
## X635               NA         NA      NA       NA
## X636               NA         NA      NA       NA
## X637               NA         NA      NA       NA
## X638               NA         NA      NA       NA
## X639               NA         NA      NA       NA
## X640               NA         NA      NA       NA
## X641               NA         NA      NA       NA
## X642               NA         NA      NA       NA
## X643               NA         NA      NA       NA
## X644               NA         NA      NA       NA
## X645               NA         NA      NA       NA
## X646               NA         NA      NA       NA
## X647               NA         NA      NA       NA
## X648               NA         NA      NA       NA
## X649               NA         NA      NA       NA
## X650               NA         NA      NA       NA
## X651               NA         NA      NA       NA
## X652               NA         NA      NA       NA
## X653               NA         NA      NA       NA
## X654               NA         NA      NA       NA
## X655               NA         NA      NA       NA
## X656               NA         NA      NA       NA
## X657               NA         NA      NA       NA
## X658               NA         NA      NA       NA
## X659               NA         NA      NA       NA
## X660               NA         NA      NA       NA
## X661               NA         NA      NA       NA
## X662               NA         NA      NA       NA
## X663               NA         NA      NA       NA
## X664               NA         NA      NA       NA
## X665               NA         NA      NA       NA
## X666               NA         NA      NA       NA
## X667               NA         NA      NA       NA
## X668               NA         NA      NA       NA
## X669               NA         NA      NA       NA
## X670               NA         NA      NA       NA
## X671               NA         NA      NA       NA
## X672               NA         NA      NA       NA
## X673               NA         NA      NA       NA
## X674               NA         NA      NA       NA
## X675               NA         NA      NA       NA
## X676               NA         NA      NA       NA
## X677               NA         NA      NA       NA
## X678               NA         NA      NA       NA
## X679               NA         NA      NA       NA
## X680               NA         NA      NA       NA
## X681               NA         NA      NA       NA
## X682               NA         NA      NA       NA
## X683               NA         NA      NA       NA
## X684               NA         NA      NA       NA
## X685               NA         NA      NA       NA
## X686               NA         NA      NA       NA
## X687               NA         NA      NA       NA
## X688               NA         NA      NA       NA
## X689               NA         NA      NA       NA
## X690               NA         NA      NA       NA
## X691               NA         NA      NA       NA
## X692               NA         NA      NA       NA
## X693               NA         NA      NA       NA
## X694               NA         NA      NA       NA
## X695               NA         NA      NA       NA
## X696               NA         NA      NA       NA
## X697               NA         NA      NA       NA
## X698               NA         NA      NA       NA
## X699               NA         NA      NA       NA
## X700               NA         NA      NA       NA
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 39 and 0 DF,  p-value: NA

A linear model doesn’t give us any significant results that we can use to say something about the predictors because it shows null p-value, null adjusted R-squared and only calculate the estimate for 39 variables + the intercept.

Extract the training and test samples

cook.app=cook[1:40,]
cook.test=cook[41:72,]

Ridge regression

Ici on a les paramètres \(\beta\) en fonction de la régularisation

library(MASS)
 plot(lm.ridge(sucrose  ~ ., data=cook.app,
   lambda = seq(0,1,0.01)))

## At lambda=0, the coefficients are large and unstable. As lambda increases, the coefficients converge toward zero

Régression ridge avec plusieurs paramètres de régularisation

 plot(lm.ridge(sucrose  ~ ., data=cook.app,
   lambda = seq(0,0.2,0.001)))

cook.ridge=lm.ridge(sucrose  ~ ., data=cook.app,
   lambda = seq(0,0.2,0.001))

Choix du paramètre de régularisation par validation croisée

The \(n=40\) observations in the training sample are randomly divided into 4 subsamples of size 10.

library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:corrplot':
## 
##     corrplot
## The following object is masked from 'package:stats':
## 
##     loadings
set.seed(87)
cvseg=cvsegments(nrow(cook.app),k=4,type="random")

Selection of the regularization parameter by cross-validation

For \(\tau\) (here \(\lambda\)) between \(0\) and \(1\), we estimate \(\beta\) and the corresponding PRESS on each subsample. We sum the 4 PRESS values to obtain the cross-validation PRESS and choose the \(\tau\) that minimizes this PRESS.

choix.kappa=function(kappamax,cvseg,nbe=100){
press=rep(0,nbe)
for (i in 1:length(cvseg)){
valid=cook.app[unlist(cvseg[i]),]
modele=lm.ridge(sucrose~.,
   data=cook.app[unlist(cvseg[-i]),],
lambda=seq(0,kappamax,length=nbe))
coeff=coef(modele)
prediction=matrix(coeff[,1],nrow(coeff),
   nrow(valid))+coeff[,-1]%*%
   t(data.matrix(valid[,-1]))
press=press+rowSums((matrix(valid[,1],
nrow(coeff),nrow(valid),byrow=T)-
   prediction)^2)
      }
kappaet=seq(0,kappamax,length=
   nbe)[which.min(press)]
return(list(kappaet=kappaet,press=press))
}

Creation of a function to plot residuals versus predicted Y values

plot.res=function(x,y,titre="title")
{
plot(x,y,col="blue",ylab="Residuals",
xlab="Predicted values",main=titre)
abline(h=0,col="green")
}

Display of results with the optimal PRESS

# execution

res=choix.kappa(0.5,cvseg,nbe=1000)
plot(seq(0,0.5,length=1000),res$press,ylab="PRESS",
xlab="lambda") # optimal value

## Base on the PRESS plot,we discover predictive performance occurs when the lambda is very close to 0. This means that the ridge penalty does not improve our prediction much in this case. The model works best with little or no regularization.

Ridge regression with the optimal parameter

kappaet=res$kappaet
cook.ridgeo=lm.ridge(sucrose~.,data=cook.app,
lambda=kappaet)
coeff=coef(cook.ridgeo)

# Computation of predicted responses

fit.rid=rep(coeff[1],nrow(cook.app))+
as.vector(coeff[-1]%*%t(data.matrix(cook.app[,-1])))

Plot of predictions versus true values

# True responses versus predicted responses

plot(fit.rid,cook.app[,"sucrose"],ylab="True responses",
xlab="Predicted responses")

# The line y = x

lines(cook.app[,"sucrose"],cook.app[,"sucrose"])

## By looking at the plot, the model predicts sucrose almost prefectly because the predicted values closely match the true responses values. This shows an excellent in-sample fit and that the chosen ridge penalty is appropriate.

Residuals

# Computation of residuals

res.rid=fit.rid-cook.app[,"sucrose"]

# Residuals versus predicted responses

plot.res(fit.rid,res.rid,titre="")

## The points are scattered randomly around the zero line. This means that the model’s predictions are generally accurate and unbiaised.

Prediction on the test sample and computation of prediction error

Since the data are centered and scaled, the appropriate transformation must be applied to obtain the estimated response.

# prediction on the test sample

ychap=rep(coeff[1],nrow(cook.test))+
as.vector(coeff[-1]%*%t(data.matrix(cook.test[,-1])))

Prediction on the test sample and computation of prediction error

plot(ychap,cook.test[,1],ylab="True responses",
xlab="Predicted responses")

## Prediction error

mean((cook.test[,1]-ychap)^2)
## [1] 9.882725

The predicted and true values on the test data follow the same trend, and this show that the model generally predict well with only moderate errors.

Comparison between Ridge and OLS coefficients

# Linear model by OLS: lambda = 0

coefflm <- coef(lm.ridge(sucrose~.,data = cook.app,
lambda=0))

Continuation

# Plot of coefficients from both models

matplot(t(rbind(coeff,coefflm)),type="l",col=c(2,1))
legend("topleft", col=c(2,1), lty=c(1,2),
legend=c("OLS beta", "Ridge beta"), cex=0.8)

## The plot shows that ridge regression smooths and stabilizes the regression coefficients compared to the OLS, helping to prevent overfitting.

LASSO: Selection of the regularization parameter

We choose 1000 values of \(\lambda\) and search, via cross-validation (on the mean squared error = MSE, using lars), for the optimal value. The plot shows the mean squared error as a function of the constraint fraction relative to \(\lambda\), \(\frac{\lambda}{|\hat{\beta}|_1}\).

# A sequence of lambda values

library(lars)
## Loaded lars 1.3
lambdas <- seq(from = 0, to = 1, length = 1000)

Continuation

# Optimal lambda or delta by CV using lars

mse.cv <- cv.lars(data.matrix(cook.app[,-1]), cook.app[,1],K = 4,se=F,index=lambdas,use.Gram=F)

lambdas.op <-lambdas[which.min(mse.cv$cv)]

From the plot, we can say the optimal lambda which minimize the cross-validation MSE is “lambda=1”.

Prediction using the optimal parameter

cook.lasso=lars(data.matrix(cook.app[,-1]),
   cook.app[,1],use.Gram=F)
fit.lasso=predict(cook.lasso,
   data.matrix(cook.app[,-1]),s=lambdas.op,
mode="fraction")$fit

Prediction using the optimal parameter

# Plot predictions

plot(fit.lasso,cook.app[,"sucrose"], ylab="True responses",
xlab="Predicted responses")

res.lasso=fit.lasso-cook.app[,"sucrose"]

Since, the points are very closed to the line y=x, we can say that the LASSO model’s predictions are very accurate.

Residuals using the optimal parameter

# Plot predictions versus residuals

plot.res(fit.lasso,res.lasso)

## The points are scattered randomly around the zero line. This means that the model’s predictions are generally accurate and unbiaised.

Prediction error using the optimal parameter

# Mean squared prediction error

pred.lasso=predict(cook.lasso,
data.matrix(cook.test[,-1]),
s=lambdas.op,mode="fraction")$fit
mean((pred.lasso-cook.test[,"sucrose"])^2)
## [1] 2.672186

Regression PCR

We choose the number of components by cross-validation.

# Mean squared errors

library(pls)
cook.pcr=pcr(sucrose~.,ncomp=28,data=cook.app,
validation="CV",segments=cvseg)
msepcv.pcr=MSEP(cook.pcr,estimate= c("train","CV"))

PCR regression

Mean errors as a function of the number of components.

plot(msepcv.pcr,type="l")

ncompo=which.min(msepcv.pcr$val["CV",,])-1
ncompo
## 6 comps 
##       6

The training MSE(in red) decreases steadily as the number of components increases because adding more components always helps fit the training data better. The plot of number of components tends to approximately 6 components, which the PCR should use to get the best predictive performance.

Regression PCR

Evolution of the proportion of variance of the variables explained by each component.

# Mean squared errors

cook.pcro=pcr(sucrose~.,ncomp=ncompo,
data=cook.app)

Regression PCR

# Mean squared errors
plot(explvar(cook.pcro),type="l",main="", pin=0.5)

## The plot shows us that the spectral data is highly redundant, which is why one component captures 83% and the others are approximately 0%.

Regression PCR

We compute the prediction error on the test sample.

fit.pcr=predict(cook.pcro,ncomp=1:6)[,,6]
res.pcr=fit.pcr-cook.app[,"sucrose"]

Prediction error on the test sample

plot.res(fit.pcr,res.pcr)

ychap=predict(cook.pcro,newdata=
cook.test)[,1,ncompo]
mean((cook.test[,"sucrose"]-ychap)^2)
## [1] 1.443031

The residual from the PCR is : “1.443031”

Regression PLS

PLS regression with up to 28 orthogonal components and residual plot with 28 components. We choose the number of components by cross-validation.

cook.pls=plsr(sucrose~.,ncomp=28,
data=cook.app,validation="CV",segments=cvseg)
print(cook.pls)
## Partial least squares regression, fitted with the kernel algorithm.
## Cross-validated using 4 random segments.
## Call:
## plsr(formula = sucrose ~ ., ncomp = 28, data = cook.app, validation = "CV",     segments = cvseg)

Regression PLS

plot(cook.pls)

msepcv.pls=MSEP(cook.pls,estimate=
   c("train","CV"))

Regression PLS

We plot the evolution of the mean squared error as a function of the number of PLS components.

plot(msepcv.pls,type="l")

## The MSE tends to decrease from 0 to 4 components and flactuate at 5 components. It stabilize around 12 components.

The PLS is more precise than the PCR because it was able to minimize the MSE with only 4 components compared to 6 components chosen by PCR.

Regression PLS

We plot the evolution of the proportion of variance of the variables explained by each component.

plot(explvar(cook.pls),type="l",main="")

## The plot shows that 2 components capture around 83% of the variance with the PLS while the PCR got only 1 component for 83%.

PLS optimal by cross-validation

## 5 comps 
##       5

##Prediction error on the test sample

PLS is better than PCR.

# mean squared error

ychap=predict(cook.plso,newdata=
cook.test)[,1,ncompo]
plot(ychap,cook.test[,1])

mean((cook.test[,"sucrose"]-ychap)^2)
## [1] 1.332782

The residual with PLS is : “1.332782” which is better than “1.443031” by the PCR.

In conclusion, PLS is better than PCR to predict quantity of sucrose in biscuits in this case.